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Analyzing coexistence and survival scenarios of Lotka-Volterra (LV) networks in which the to- 
tal biomass is conserved is of vital importance for the characterization of long-term dynamics of 
ecological communities. Here, we introduce a classification scheme for coexistence scenarios in 
these conservative LV models and quantify the extinction process by employing the Pfafhan of the 
network's interaction matrix. We illustrate our findings on global stability properties for general 
systems of four and five species and find a generalized scaling law for the extinction time. 
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Understanding the stability of ecological networks is of 
pivotal importance in theoretical biology [1, 2]. Coexis- 
tence and extinction of species depend on many factors 
such as inter- and intra-species interactions [3, 4], pop- 
ulation size [5-9], and mobility of individuals [10-16]. 
An intriguing question is how the stability of ecosystems 
depends on the interaction network between species. Is 
it the topology of the network (whose links may arise 
through predation, competition over common resources, 
or mutual cooperation) that sets the level of biodiversity? 
And how important is the strength of a single interaction 
link? Stable coexistence can, for example, be observed 
for natural populations in non-hierarchical networks that 
are comprised of species that interact in a competitive 
and predator-prey like manner [17, 18]. By understand- 
ing the interplay between the structure of the interaction 
network and the strengths of its links, it is possible to 
reveal mechanisms that underlie this stability. 

A paradigm in addressing these ecologically impor- 
tant questions from a theoretical perspective are Lotka- 
Volterra (LV) models [19, 20] in which the total biomass 
of species is conserved. These conservative LV sys- 
tems [12, 21, 22] originate in the well-mixed limit from 
agent-based formulations of reaction-diffusion systems, 
where individuals of S different species A\, A2, ■ . . ,As 
compete directly with each other following the simpli- 



— > Ai + Ai. Species 
and immediately re- 



fied reaction scheme [23] : A t + Aj 
Ai beats species Aj with rate ki 3 
places an individual of species Aj with an own offspring. 
Species Aj is thus degraded at the same rate such that 
the interaction matrix Gs = {kij}i,j is skew-symmetric. 
The interaction network can be visualized by a graph; see 
Fig. 1. Neglecting demographic fluctuations [24], the de- 
terministic dynamics for the species' concentration vector 
x = (xi, . . . , xs) T is given by the rate equations (REs): 



(Gs~x)i , for alH = 1, . . . , S 



(1) 



This conservative LV model has been investigated as a 
prototype to understand principles of biodiversity from 
a theoretical point of view [8, 25]. While these sys- 
tems are also of central importance to many other fields 



of science (e.g., plasma physics [26], evolutionary game 
theory [27, 28], and chemical kinetics [29]), no general 
scheme to classify coexistence, survival, and extinction 
of species has been established so far. It is frequently as- 
sumed that the topology of the interaction network alone 
determines coexistence of species [30, 31], i.e., that such 
systems can be regarded as Boolean networks [32]. Re- 
cent investigations of specific topologies indicate, how- 
ever, that knowledge about the network topology may 
not suffice to conclude whether all species coexist or if 
some of them go extinct [33-35]. These questions on 
global stability properties have been previously addressed 
successfully for various particular LV systems [27, 36] and 
for hierarchical networks [37, 38]. 

In this letter, we present a general classification of co- 
existence scenarios in conservative LV networks with an 
arbitrary number of species. We elucidate the conse- 
quences of the interplay between the network structure 
and the strengths of its interaction links on global sta- 
bility. By analyzing conserved quantities, we find condi- 
tions on the reaction rates that yield coexistence of all 
species. In our mathematical framework this amounts to 
the characterization of positive kernel elements of the in- 
teraction matrix: By employing the algebraic concept of 
the Pfaffian of a skew-symmetric matrix, we are able to 
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FIG. 1. (Color online) Two interaction topologies specify- 
ing the conservative LV systems, (a) The general cyclic four 
species systems (4SS). (b) The general cyclic five species sys- 
tem (5SS) as a natural extension of the rock-papers-scissors 
configuration (RPS) [39]. 
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generalize previous approaches [34, 40] and to quantify 
the extinction process when no conserved quantities ex- 
ist. We illustrate our general results for coexistence and 
survival scenarios of four and five species systems (4SS 
and 5SS), cf. Fig. 1. Moreover, we demonstrate the im- 
plications of our findings for the stability of stochastic 
systems: We show how the extinction time diverges with 
the distance to the critical rate at which coexistence of 
all species is observed. 

First, we discuss some general results for the REs (1) 
before the specific interaction topologies in Fig. 1 are ana- 
lyzed. In order to characterize the stability of the generic 
LV system, we study conserved quantities. We elabo- 
rate on the form of conserved quantities, under which 
conditions they exist at all, and how many conserved 
quantities there are for a given interaction network. 
Since the interaction matrix Gs is skew-symmetric, the 
REs (1) conserve the sum over all species' concentrations 
t = xi + . . . Xs, independent of the interaction scheme. 
Hence, the dynamics can be normalized onto the (S — 1)- 
dimensional simplex where all concentrations are non- 
negative and add up to 1. The vertices of the simplex 
correspond to the extinction of all but one species, its 
edges reflect the extinction of all but two species, and 
so on. Further conserved quantities have previously been 
derived [20, 27, 40]. Interestingly, these 

conserved quantities can be obtained from solutions of 
the linear problem GsP — because f = — r(Gsp,x), 
with p = (pi, . . . ,ps) T ■ One infers that r is conserved if 
the exponent vector p is an eigenvector corresponding to 
eigenvalue [40], or in other words, if p lies in the kernel 
of the matrix Gs- 

Coexistence means that all concentrations stay away 
from the boundary of the simplex by a finite distance for 
all times. Since the species' concentrations are bounded 
to the interval [0, 1], one concludes from the structure of 
the conserved quantity r that all S species coexist if the 
kernel of the interaction matrix is positive, i.e., one finds 
an element p in the kernel of Gs with positive entries 
Pi > for all i. Hence, to reveal coexistence scenarios 
in the conservative LV model, one has to characterize 
the kernel of the interaction matrix Gs and identify its 
positive elements. Note that this conclusion goes beyond 
stating that a positive kernel element corresponds to a 
stationary point in the inside of the simplex; see REs (1). 

The existence of conserved quantities constrains the 
dynamics to a submanifold of the simplex whose dimen- 
sion D c is determined as follows. The rank of a skew- 
symmetric matrix is always even, because its non-zero 
eigenvalues are purely imaginary, conjugate pairs. The 
rank-nullity theorem [41] then implies that the dimension 
of the kernel of Gs is odd whenever S is odd, and even 
whenever S is even. Each linearly independent kernel 
element gives rise to an independent conserved quantity 
r which constrains the degrees of freedom of the tra- 
jectory. Together with tq, one finds that the dynamics 
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FIG. 2. (Color online) Coexistence and survival in the general 
cyclic 4SS are controlled by the Pfaffian of the interaction 
matrix, (a) For Pf (G4) = 0, one obtains coexistence of all 
species on periodic orbits, (b) Deterministic survival diagram: 
for Pf (G4) < 0, species A, B, and D survive in a stable RPS 
configuration, whereas A, G, and D survive for Pf (G4) > 0. 



in case of non-stationary motion is constrained to a de- 
formed sphere of dimension D c = S — 1 — dimKer (Gs) 
for a positive kernel; see the Supplemental Material 
(SM) for mathematical details. Thus, coexistence in 
high-dimensional systems is generically observed on non- 
periodic trajectories (D c > 1); see Movie Ml of SM. Only 
if the reaction rates are fine-tuned to a positive and max- 
imal kernel of dimension 5 — 2, the dynamics is restricted 
to periodic orbits (D c = 1); see Fig. 2(a) and Movie M2 
of SM. In particular, for S = 3 or 4, a positive kernel 
immediately implies coexistence on periodic orbits. This 
follows from the fact that with three species, the kernel is 
always one-dimensional. For the general 4SS, the dimen- 
sion of the kernel of the interaction matrix is either or 2. 
A two-dimensional, positive kernel yields coexistence on 
periodic orbits; see Fig. 2(a). If dimKer (Gs) = 0, i.e., if 
the kernel is trivial, one observes extinction of species as 
detailed below. 

Next, we focus on the mapping between the reaction 
rates in Gs and its kernel elements in order to find the 
stationary points. To this end, we apply the concepts of 
the Pfaffian and of the adjugate matrix [41, 42]. The 
Pfaffian is a simpler form of the determinant tailored 
to skew-symmetric matrices with the property that its 
square equals the value of the determinant. In contrast 
to the non-negative determinant of skew-symmetric ma- 
trices, the Pfaffian carries a sign which will turn out to be 
crucial for our purposes. For a skew-symmetric matrix, 
the Pfaffian can be computed recursively as: 



Pf(G s ) =£(-!)'. fc M .Pf(G B ) 



(2) 



1=2 



with Gn being the matrix where both the first and i- 
th column and row have been removed from the matrix 
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Gs- The Pfaffian of a 2 x 2 skew-symmetric matrix G 2 — 
{kAB}, is given by Pf (G 2 ) = ^ab- For the interaction 
matrix corresponding to the LV network in Fig. 1(a), 

f k AB k AC -k DA \ 

q ~k A B k B c k B D 

-k AC -k BC k C D 

\ k DA -k BD -k C o / 

the PfafBan is Pf (G 4 ) = k AB k C D - k AC k BD - k DA k B c- 

The Pfaffian always vanishes for odd S as opposed to 
systems with an even number of species [42] . In the latter 
case, the Pfaffian is zero only if a constraint on the reac- 
tion rates is fulfilled. If the Pfaffian vanishes, one finds 
more kernel elements than just the null vector and, thus, 
conserved quantities of form r exist. In the following, we 
distinguish between even and odd S. 

For an even number of species and a two-dimensional 
kernel, positive kernel elements can be identified via the 
adjugate matrix Rs which is a generalized inverse of the 
interaction matrix such that Gs ■ Rs = — Pf(Gs) • Is, 
with I5 being the identity matrix [42]. The adjugate 
matrix can be computed as (Rs)ij = (— l) CT Pf (Gy) 
where (— l) a denotes the sign of the permutation a = 
(i j 1 . . . i . . . j . . . S) , and the columns of Rs give two 
independent kernel elements of Gs- 

As an example, consider again the general cyclic 4SS 
depicted in Fig. 1(a). By setting all reaction rates equal 
to each other (e.g., to 1), the Pfaffian does not vanish 
and, therefore, not all species can coexist. Only when the 
rates are chosen such that Pf (G4) — 0, do we obtain two 
independent kernel elements of G 4 : From its adjugate 
matrix, i? 4 , we identify pi = (kcD, 0, k BAl k A c) T and 
P2 = (k B D,k B A,0,k AB ) T . We infer the two conserved 
quantities n = x A CD x k ( f A x D AC and r 2 = x A BD x k B DA x k D AB , 
and conclude that the kernel is positive and coexistence 
occurs on periodic orbits; see Fig. 2(a). Hence, classifying 
LV networks in terms of their topology is incomplete; the 
strengths of the interaction links are crucial in general. 

In general, if the Pfaffian for a system with even S is 
non-zero, i.e., when only the null vector lies in the ker- 
nel, coexistence of all species is not possible. Still one 
can quantify the extinction process by generalizing an 
approach of Durney et al. [34] for a system with S = 4 to 
systems composed of an arbitrary even number of species. 
We define the function p — x' 1 . . . x q s s in the same way 
as the conserved quantity t, but this time choosing the 
exponent vector = — R$l with 1 = (1, . . . , 1) T . It is 
straightforward to show that this function evolves expo- 
nentially in time: 

p(t) = p(0) ■ e - pf (^H , (3) 

generalizing previous investigations [24, 33-35]. It is 
quite remarkable that p quantifies the global collective 
dynamics of systems with an arbitrary interaction topol- 
ogy and even S. Depending on the sign of the Pfaffian, 



p grows or decays exponentially fast with the Pfaffian of 
the interaction matrix as rate. Since the system's dy- 
namics is driven towards the boundary of the simplex, 
one can conclude on the extinction of some species. This 
feature of p is reminiscent of a Lyapunov function; note 
also that p becomes a conserved quantity r if the Pfaffian 
is zero. An interesting question for future investigations 
is to ask whether further quantities exist that character- 
ize the dynamics of conservative LV networks. 

For the general 4SS shown in Fig. 1(a), we find q 4 = 
(-k CD + k BD - k BC , k CD + k DA + k AC , -k BD - k DA - 
k AB , k BC — k AC + k AB ) T . The fact that (q 4 ) 2 is always 
positive suggests that species B goes extinct for a posi- 
tive Pfaffian, and that the converse holds true for (q 4 )3 
and species G for a negative Pfaffian. In both cases, the 
system tends to a stable rock-paper-scissors (RPS) con- 
figuration. In summary, we derive the survival diagram 
shown in Fig. 2(b). Interestingly, A and D always sur- 
vive in this topology although D can be easily tuned to 
be the weakest species. We emphasize that this result de- 
pends on the sign of the Pfaffian and cannot be obtained 
from applying the concept of the determinant. Again, 
since the Pfaffian of the interaction matrix characterizes 
the dynamics of this 4SS, its topology alone does not 
determine the long-time dynamics. These findings unify 
previous results for other 4SS [24, 33, 34], and show that 
rules like "survival of the strongest" or "survival of the 
weakest" [25, 43] cannot be formulated in general. 

For an odd number of species, the kernel of Gg is al- 
ways nontrivial. In general, if dimkerGg = 1, we deter- 
mine the independent kernel element via the adjugate 
vector [42], r s = (Pf (G t ) , -Pf (G 2 ) , . . . , Pf (G S )) T , 
which enables us to investigate the influence of the re- 
action rates on the survival scenarios. For 5 = 3, only 
the well-studied RPS topology [8, 27] leads to a positive 
adjugate vector r 3 . In other words, coexistence of all 
three species depends only on the topology of the net- 
work. This behavior is unique to S = 3 and changes 
dramatically for systems with more than three species. 

We illustrate the importance of the reaction rates for 
a system of five interacting species; see Fig. 1(b). This 
interaction topology where each species dominates two 
species and is outperformed by the two remaining species, 
recently gained attention as a natural extension of the 
RPS game [30, 39, 44]. For specificity, we investigate 
the dependence of the survival scenarios on the rate k AB 
with which species A beats species B and chose the other 
rates (see Fig. 3(b), left inset) such that either five or 
four species survive depending on the value of k AB ; sec 
Fig. 3(a). The kernel of the interaction matrix depends 
on k AB and is characterized by the adjugate vector r 5 = 
(0,0,3/cab - 15,5 - k AB , 5k AB - 25) T . For k AB ± 5, 
the kernel is one-dimensional and non-positive, and four 
species survive. In contrast, for k AB — 5, r 5 equals the 
null vector which in turn means that the kernel becomes 
three-dimensional [42]. Since we have ensured that the 
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FIG. 3. (Color online) Stability of the cyclic 5SS. (a) For the 
interaction scheme (left inset of (b)), one obtains coexistence 
of all species for the critical rate kAB = 5. (b) Stability of 
the stochastic system, reflected by the extinction time T ext , 
peaks at the critical rate, which becomes more pronounced 
as TV — > oo. We find a scaling law for T ext in the distance to 
the critical rate (right inset). Initial conditions were chosen 
as x(0) = 1/5-1. Larger line gap corresponds to smaller N. 



kernel is also positive, we obtain coexistence of all five 
species on periodic orbits (D c = 1). 

Finally, we discuss the implications of our findings by 
asking how demographic noise affects the stability of 
stochastic LV systems. We analyze ecological LV sys- 
tems with a finite number N of interacting individuals 
in the eye of the knowledge gained from the determinis- 
tic analysis. It has been shown that due to demographic 
fluctuations the system ultimately reaches an absorbing 
state that is characterized by the extinction of all but 
one species [45-48]. Moreover, the scaling behavior of 
the mean extinction time with the system size N charac- 
terizes the stability of the interaction network [14, 47]. 

As an example, we continue the discussion of the 5SS 
from Fig. 3(b), left inset. We have carried out exten- 
sive computer simulations employing the Gillespie algo- 
rithm [49] to measure the time T cxt until the first species 
has become extinct for different system sizes N and dif- 
ferent reaction rates kAB- The results are displayed in 
Fig. 3(b) and highlight the significance of the determinis- 
tic drift underlying the stochastic dynamics. We observe 
a peak in the extinction time as the reaction rate kAB 
approaches the critical value fc cr = 5 for which we obtain 
coexistence of all species in the deterministic case. The 
divergence of the extinction time for kAB ~> k cr becomes 
more pronounced for larger system sizes as the system 
reaches the deterministic limit for N — > oo. 

A scaling analysis reveals how the extinction time 
peaks in the vicinity of the coexistence scenario. Near 
the critical rate, the extinction time scales linearly with 
the system size leading to neutrally stable interaction 



networks [8, 24, 50]. At larger distance from the criti- 
cal rate, the deterministic driving force to the absorbing 
boundary becomes more dominant than the demographic 
fluctuations; see Fig. 3(b), right inset. The interplay be- 
tween the stochastic and deterministic forces is reflected 
by the scaling law: 



N 



In N 



-fccrl 



for k A B = k c 
for k A B k c 



(4) 



which extends the linear scaling T cxt oc N of neutral co- 
existence. We observe a power-law dependence in the 
distance of the reaction rates to the critical rate and log- 
arithmic scaling with N for attracting boundaries [8, 51]. 

The observed scaling law (4) for kAB ^ k CI can be 
attributed to the exponentially fast extinction of species 
Xi = Xi(0) exp (—ait); see Eq. (1). The extinction rate a.i 
is computed via the temporal average over the trajectory 
(x) as ctj = — (Gs (x))i, which becomes linear in the 
distance to the critical rate \kAB — k CI \ for large times. 
The logarithmic dependence on N follows by defining 
that a species with concentration xi less than 1/N has 
become extinct. With this scaling behavior at hand, we 
are able to compare the ecological stability of different 
interaction networks based on our analysis of the REs (1). 

In this Letter, we investigated global stability proper- 
ties of conservative LV networks. By employing the Pfaf- 
fian of the interaction matrix, we revealed the relation 
between the reaction rates and the conditions for coexis- 
tence, and exemplified the implications for the stability 
of ecological networks with finite populations. We expect 
that our results will also stimulate further progress for the 
investigation of extinction scenarios. Beyond analyzing 
whether an ecosystem is stable or unstable, it would be 
highly interesting to actually predict which of its species 
ultimately survive for a general conservative LV system. 
This would, for example, allow us to predict the eventual 
outcome of an unstable version of the five species system 
shown in Fig. 1(b), and to formulate the conditions under 
which 3- or 4-species cycles are attained. First insights 
into these extinction dynamics will be outlined in a future 
publication [52]. We believe that a full characterization 
of general conservative LV dynamics is possible. 
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